import sympy as sp

# S1 -> S2 
# 在S1中的速度为 (b1x,0,0)
# b = v/c
g1,b1x = sp.symbols('g1,b1x',positive=True);

# S1 -> S3
# 在S1中的速度为 (b1x+dbx, dby,0)
dbx,dby = sp.symbols('dbx,dby',positive=True);

# S1 -> S2 Lorentz变换矩阵
L2 = sp.Matrix([[g1,-g1*b1x,0,0],[-g1*b1x,g1,0,0],[0,0,1,0],[0,0,0,1]])
L2_inv = sp.Matrix([[g1,g1*b1x,0,0],[g1*b1x,g1,0,0],[0,0,1,0],[0,0,0,1]])


# S1 -> S3 Lorentz变换矩阵
# 根据任意方向的Lorentz变换矩阵公式（i,j=1,2,3） L=[g,-g*bi;-g*bi,delta_{ij} + (g-1)*bi*bj/(\sum bi^2)]
# 取一阶小量
L3_00 = g1 + g1**3*b1x*dbx
L3_01 = -(g1*b1x+g1**3*dbx)
L3_02 = -(g1*dby)
L3_11 = L3_00
L3_12 = (g1-1)/b1x*dby

L3 = sp.Matrix(
[   
    [L3_00,L3_01,L3_02,0],
    [L3_01,L3_11,L3_12,0],
    [L3_02,L3_12,1,    0],
    [0    ,0    ,0,    1],
]
)
print('L2:')
print(sp.latex(L2))
print('L3:')
print(sp.latex(L3))

# S2->S3 变换矩阵（非Lorentz变换矩阵）
B = L3@L2_inv

print('B:')
print(sp.latex(B))

print('B_asym')
B_asym = B-B.T
B_asym = B_asym.applyfunc(sp.simplify)
print(B_asym)

# B = L R
# 应用旋转，得到L
dTheta = sp.symbols('dTheta',positive=True)
R = sp.Matrix([[1,0,0,0],[0,1,dTheta,0],[0,-dTheta,1,0],[0,0,0,1]])
L = B@R.T

print('L:')
print(sp.latex(L))

print('dTheta:')
expr = L[2,1]-L[1,2]
expr = expr.expand()
expr = -(expr-expr.coeff(dTheta)*dTheta)/expr.coeff(dTheta)
expr = expr.simplify()
print(expr)
